/*
 *    Geotoolkit.org - An Open Source Java GIS Toolkit
 *    http://www.geotoolkit.org
 *
 *    (C) 1999-2012, Open Source Geospatial Foundation (OSGeo)
 *    (C) 2009-2012, Geomatys
 *
 *    This library is free software; you can redistribute it and/or
 *    modify it under the terms of the GNU Lesser General Public
 *    License as published by the Free Software Foundation;
 *    version 2.1 of the License.
 *
 *    This library is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *    Lesser General Public License for more details.
 *
 *    NOTE: permission has been given to the JScience project (http://www.jscience.org)
 *          to distribute this file under BSD-like license.
 */
package org.geotoolkit.nature;

import static java.lang.Math.*;


/**
 * Sea water properties as a function of salinity, temperature and pressure.
 * Density is computed using the 1980 definition of Equation of State (EOS80).
 * Units are:
 * <p>
 * <ul>
 *   <li>Salinity: Pratical Salinity Scale 1978 (PSS-78).</li>
 *   <li>Temperature: Celsius degrees according International Temperature Scale 1968 (ITS-68).</li>
 *   <li>Pressure: decibars (1 dbar = 10 kPa).
 * </ul>
 *
 * @author Bernard Pelchat (MPO)
 * @author Martin Desruisseaux (MPO, IRD)
 * @version 3.00
 *
 * @since 1.0
 * @module
 */
public final class SeaWater {
    /*
     * Note: Les algorithmes originaux de l'UNESCO recevaient en entrés
     *       des pressions en décibars. Les algorithmes écrites par
     *       Bernard Pelchat recevaient en entrés des pressions en
     *       MegaPascal. La première ligne de code des algorithmes
     *       de Bernard Pelchat multipliait donc les pressions par
     *       100, afin de les convertir en decibars.
     */

    /**
     * Conductivity (in mS/cm) of a standard sea water sample.
     * S is for <cite>Siemens</cite> (or Mho, its the same...).
     */
    public static final double STANDARD_CONDUCTIVITY = 42.914;

    /**
     * Coéfficients de l'équation d'état EOS-80. La densité
     * calculée par ces coéfficients est la densité Sigma-T.
     */
    private static final double
        EOS80_A[]  = {-28.263737,    6.793952E-2, -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9},
        EOS80_B[]  = {  8.24493E-1, -4.0899E-3,    7.6438E-5,  -8.2467E-7,     5.3875E-9},
        EOS80_C[]  = { -5.72466E-3,  1.0227E-4,   -1.6546E-6},
        EOS80_D    =    4.8314E-4,
        EOS80_E[]  = {-1930.06,      148.4206,    -2.327105,    1.360477E-2, -5.155288E-5},
        EOS80_F[]  = {   54.6746,     -6.03459E-1, 1.09987E-2, -6.1670E-5},
        EOS80_G[]  = {    7.944E-2,    1.6483E-2, -5.3009E-4},
        EOS80_H[]  = {   -1.194975E-1, 1.43713E-3, 1.16092E-4, -5.77905E-7},
        EOS80_I[]  = {    2.2838E-3,  -1.0981E-5, -1.6078E-6},
        EOS80_J    =      1.91075E-4,
        EOS80_K[]  = {    3.47718E-5, -6.12293E-6, 5.2787E-8},
        EOS80_M[]  = {   -9.9348E-7,   2.0816E-8,  9.1697E-10},
        EOS80_N[]  = {21582.27,        3.359406,   5.03217E-5},
        RHO_35_0_0 =  1028.1063,
        DR_35_0_0  =  28.106331;

    /**
     * Coéfficients de l'équation d'état EOS-80. La densité
     * calculée par ces coéfficients est la densité "vrai".
     */
    private static final double
        EOS80_At[] = {999.842594,   6.793952E-2, -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9},
        EOS80_Et[] = {19652.21,   148.4206,      -2.327105,    1.360477E-2, -5.155288E-5},
        EOS80_Ht[] = { 3.239908,    1.43713E-3,   1.16092E-4, -5.77905E-7},
        EOS80_Kt[] = { 8.50935E-5, -6.12293E-6,   5.2787E-8};

    /**
     * Coéfficients de l'équation de la salinité PSS-78.
     */
    private static final double
        PSS78_A[] = { 0.0080,    -0.1692,    25.3851,      14.0941,   -7.0261,  2.7081},
        PSS78_B[] = { 0.0005,    -0.0056,    -0.0066,      -0.0375,    0.0636, -0.0144},
        PSS78_C[] = { 0.6766097,  2.00564E-2, 1.104259E-4, -6.9698E-7, 1.0031E-9},
        PSS78_D[] = { 3.426E-2,   4.464E-4,   4.215E-1,    -3.107E-3},
        PSS78_E[] = { 2.070E-5,  -6.370E-10,  3.989E-15},
        PSS78_G[] = {-0.1692,    50.7702,    42.2823,     -28.1044, 13.5405},
        PSS78_H[] = {-0.0056,    -0.0132,    -0.1125,       0.2544, -0.0720},
        PSS78_K   =   0.0162;

    /**
     * Coéfficients pour les salinités élevées,
     */
    private static final double[]
        PSS78_AR = {  7.737,   -9.819,     8.663, -2.625},
        PSS78_AT = {  3.473E-2, 3.188E-3, -4.655E-5},
        PSS78_CR = {-10.01E-2,  4.82E-2,  -6.682E-4};

    /**
     * Constantes nécessaires au calcul de la chaleur spécifique.
     *
     * @see #specificHeat
     */
    private static final double
        HEAT_AA[] = {-7.643575,   0.1072763,  -1.38385E-3},
        HEAT_BB[] = { 0.1770383, -4.07718E-3,  5.148E-5},
        HEAT_CC[] = { 4217.4,    -3.720283,    0.1412855,  -2.654387E-3, 2.093236E-5},
        HEAT_A[]  = {-4.9592E-1,  1.45747E-2, -3.13885E-4,  2.0357E-6,   1.7168E-8},
        HEAT_B[]  = { 2.4931E-4, -1.08645E-5,  2.87533E-7, -4.0027E-9,   2.2956E-11},
        HEAT_C[]  = {-5.422E-8,   2.6380E-9,  -6.5637E-11,  6.136E-13},
        HEAT_D[]  = { 4.9247E-3, -1.28315E-4,  9.802E-7,    2.5941E-8,  -2.9179E-10},
        HEAT_E[]  = {-1.2331E-4, -1.517E-6,    3.122E-8},
        HEAT_F[]  = {-2.9558E-6,  1.17054E-7, -2.3905E-9,   1.8448E-11},
        HEAT_G    =   9.971E-8,
        HEAT_H[]  = { 5.540E-10, -1.7682E-11,  3.513E-13},
        HEAT_J    =  -1.4300E-12;

    /**
     * Constantes nécessaires au calcul de la température adiabétique.
     *
     * @see #adiabeticTemperatureGradient
     */
    private static final double
        GRAD_A[] = { 3.5803E-05,  8.5258E-06, -6.8360E-08,  6.6228E-10},
        GRAD_B[] = { 1.8932E-06, -4.2393E-08},
        GRAD_C[] = { 1.8741E-08, -6.7795E-10,  8.7330E-12, -5.4481E-14},
        GRAD_D[] = {-1.1351E-10,  2.7759E-12},
        GRAD_E[] = {-4.6206E-13,  1.8676E-14, -2.1687E-16};

    /**
     * Constantes nécessaires au calcul de la profondeur.
     *
     * @see #depth
     */
    private static final double
        DEPTH_C[] = {9.72659, -2.2512E-5, 2.279E-10, -1.82E-15};

    /**
     * Constantes nécessaires au calcul de la vitesse du son.
     *
     * @see #soundVelocity
     */
    private static final double
        SOUND_A0[] = { 1.389,     -1.262E-2,    7.164E-5,   2.006E-6,  -3.21E-8},
        SOUND_A1[] = { 9.4742E-5, -1.2580E-5,  -6.4885E-8,  1.0507E-8, -2.0122E-10},
        SOUND_A2[] = {-3.9064E-7,  9.1041E-9,  -1.6002E-10, 7.988E-12},
        SOUND_A3[] = { 1.100E-10,  6.649E-12,  -3.389E-13},
        SOUND_B0[] = {-1.922E-2,  -4.42E-5},
        SOUND_B1[] = { 7.3637E-5,  1.7945E-7},
        SOUND_C0[] = {1402.388,    5.03711,    -5.80852E-2, 3.3420E-4,  -1.47800E-6, 3.1464E-9},
        SOUND_C1[] = {0.153563,    6.8982E-4,  -8.1788E-6,  1.3621E-7,  -6.1185E-10},
        SOUND_C2[] = {3.1260E-5,  -1.7107E-6,   2.5974E-8, -2.5335E-10,  1.0405E-12},
        SOUND_C3[] = {-9.7729E-9,  3.8504E-10, -2.3643E-12 },
        SOUND_D0   =   1.727E-3,
        SOUND_D1   =  -7.9836E-6;

    /**
     * Constantes nécessaires au calcul de la saturation en oxygène dissous.
     *
     * @see #saturationO2
     */
    private static final double
        O2_AT[] = {-135.29996, 1.572288E+5, -6.637149E+7, 1.243678E+10, -8.621061E+11},
        O2_AS[] = {0.020573, -12.142, 2363,1};




    /**
     * Do not allow instantiation of this class.
     */
    private SeaWater(){
    }

    /**
     * Computes density as a function of salinity, temperature and pressure.
     *
     * @param S Salinity PSS-78 (0 to 42)
     * @param T Temperature ITS-68 (-2 to 40°C)
     * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
     * @return  Density (kg/m³).
     */
    public static double density(final double S, final double T, double P) {
        P /= 10.0;

        // Pure water density at atmospheric pressure
        final double RHO_0_T_0 = poly(T,EOS80_At);

        // Sea water density at atmospheric pressure
        final double SR = sqrt(S);
        final double RHO_S_T_0 = (EOS80_D*S + poly(T,EOS80_C)*SR + poly(T,EOS80_B))*S + RHO_0_T_0;

        // Compression terms
        final double K_S_T_0 = (poly(T,EOS80_F) + poly(T,EOS80_G)*SR)*S + poly(T,EOS80_Et);
        final double K_S_T_P = K_S_T_0 + ((EOS80_J*SR + poly(T,EOS80_I)) * S + poly(T,EOS80_Ht) +
                               (poly(T,EOS80_Kt) + poly(T,EOS80_M) * S) * P) * P;
        return RHO_S_T_0/( 1.0 - P/K_S_T_P );
    }

    /**
     * Computes density sigma-T as a function of salinity, temperature and pressure.
     * Density Sigma-T is equivalent to the true density minus 1000&nbsp;kg/m³, and
     * has typical values around 35. This computation avoid some rouding errors
     * occuring in the true density computation.
     *
     * @param S Salinity PSS-78 (0 to 42)
     * @param T Temperature ITS-68 (-2 to 40°C)
     * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
     * @return  Density Sigma-T (kg/m³).
     */
    public static double densitySigmaT(final double S, final double T, double P) {
        P /= 10.0;
        // Sea water density at atmospheric pressure
        final double SR = sqrt(S);
        final double RHO = (EOS80_D*S + poly(T,EOS80_C)*SR + poly(T,EOS80_B))*S + poly(T,EOS80_A);

        // Specific volume at atmospheric pressure
        final double V_35_0_0    = 1.0/RHO_35_0_0;
        final double SVAN_S_T_0  = -RHO*V_35_0_0/(RHO+RHO_35_0_0);
        if (P <= 0) {
            return RHO + DR_35_0_0;
        }
        // Compression terms, DK = K(S,T,P) - K(35,0,P)
        final double K0 = (poly(T,EOS80_F) + poly(T,EOS80_G)*SR)*S + poly(T,EOS80_E);
        final double DK = K0 + (((EOS80_J * SR + poly(T,EOS80_I)) * S + poly(T,EOS80_H)) +
                          (poly(T,EOS80_K) + poly(T,EOS80_M) * S) * P) * P;

        final double K_35_0_P = poly(P,EOS80_N);
        final double V_S_T_0  = SVAN_S_T_0 + V_35_0_0;
        final double SVANS    = SVAN_S_T_0 * (1.0 - P/K_35_0_P) + V_S_T_0 * P * DK /
                                (K_35_0_P * (K_35_0_P + DK));

        // Compute density anomaly
        final double V_35_0_P  = V_35_0_0*( 1.0 - P/K_35_0_P );
        final double DR_35_0_P = P/(K_35_0_P*V_35_0_P);
        final double DVAN      = SVANS/( V_35_0_P*( V_35_0_P + SVANS ) );
        return DR_35_0_0 + DR_35_0_P - DVAN;
    }

    /**
     * Computes volume as a function of salinity, temperature and pressure.
     * This quantity if the inverse of density. This method is equivalent
     * to <code>1/{@link #density density}(S,T,P)</code>.
     *
     * @param S Salinity PSS-78 (0 to 42)
     * @param T Temperature ITS-68 (-2 to 40°C)
     * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
     * @return  Volume (m³/kg).
     */
    public static double volume(final double S, final double T, double P) {
        P /= 10.0;
        // Sea water density at atmospheric pressure
        final double SR = sqrt(S);
        final double RHO = (EOS80_D*S + poly(T,EOS80_C)*SR + poly(T,EOS80_B))*S + poly(T,EOS80_A);

        // Specific volume at atmospheric pressure
        final double V_35_0_0    = 1.0/RHO_35_0_0;
        final double SVAN_S_T_0  = -RHO*V_35_0_0/(RHO+RHO_35_0_0);
        if (P <= 0) {
            return SVAN_S_T_0 + V_35_0_0;
        }
        // Compression terms, DK = K(S,T,P) - K(35,0,P)
        final double K0 = (poly(T,EOS80_F) + poly(T,EOS80_G) * SR) * S + poly(T,EOS80_E);
        final double DK = K0 + (((EOS80_J * SR + poly(T,EOS80_I)) * S + poly(T,EOS80_H)) +
                          (poly(T,EOS80_K) + poly(T,EOS80_M) * S) * P) * P;

        final double K_35_0_P = poly(P,EOS80_N);
        final double V_S_T_0  = SVAN_S_T_0 + V_35_0_0;
        return (SVAN_S_T_0 + V_35_0_0) * (1.0 - P/K_35_0_P) + V_S_T_0 * P * DK / (K_35_0_P * (K_35_0_P + DK));
    }

    /**
     * Computes volumic anomaly as a function of salinity, temperature and pressure.
     * Volumic anomaly is defined as the sea water sample's volume minus a standard
     * sample's volume, where the standard sample is a sample of salinity 35, temperature
     * 0°C and the same pressure. In pseudo-code, {@code volumeAnomaly} is equivalent
     * to <code>{@link #volume volume}(S,T,P)-{@link #volume volume}(35,0,P)</code>.
     *
     * @param S Salinity PSS-78 (0 to 42)
     * @param T Temperature ITS-68 (-2 to 40°C)
     * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
     * @return  Volumic anomaly (m³/kg).
     */
    public static double volumeAnomaly(final double S, final double T, double P) {
        P /= 10.0;
        // Sea water density at atmospheric pressure
        final double SR = sqrt(S);
        final double RHO = (EOS80_D*S + poly(T,EOS80_C)*SR + poly(T,EOS80_B))*S + poly(T,EOS80_A);

        // Specific volume at atmospheric pressure
        final double V_35_0_0    = 1.0/RHO_35_0_0;
        final double SVAN_S_T_0  = -RHO*V_35_0_0/(RHO+RHO_35_0_0);
        if (P <= 0) {
            return SVAN_S_T_0;
        }
        // Compression terms, DK = K(S,T,P) - K(35,0,P)
        final double K0 = (poly(T,EOS80_F) + poly(T,EOS80_G)*SR)*S + poly(T,EOS80_E);
        final double DK = K0 + (((EOS80_J * SR + poly(T,EOS80_I)) * S + poly(T,EOS80_H)) +
                          (poly(T,EOS80_K) + poly(T,EOS80_M) * S) * P) * P;

        final double K_35_0_P = poly(P,EOS80_N);
        final double V_S_T_0  = SVAN_S_T_0 + V_35_0_0;
        return (SVAN_S_T_0*(1.0 - P/K_35_0_P) + V_S_T_0 * P * DK / (K_35_0_P * (K_35_0_P + DK)));
    }

    /**
     * Practical salinity scale 1978 definition
     * with temperature correction, XR = SQRT( Rt )
     */
    private static double sal(double RT, double XT) {
        return poly(RT,PSS78_A) + (XT/(1.0+PSS78_K*XT)) * poly(RT,PSS78_B);
    }

    /**
     * {@code dsal(RT,XT)} function for derivative
     * of {@code sal(RT,XT)} with <var>RT</var>.
     */
    private static double dsal(double RT, double XT) {
        return poly(RT,PSS78_G) + (XT/(1.0+PSS78_K*XT)) * poly(RT,PSS78_H);
    }

    /**
     * Computes salinity as a function of conductivity, temperature and pressure.
     *
     * @param C Conductivity in mS/cm (millisiemens by centimeters). Multiply
     *          par {@link #STANDARD_CONDUCTIVITY} if {@code C} is not a
     *          real conductivity, but instead the ratio between the sample's
     *          conductivity and the standard sample's conductivity.
     * @param T Temperature ITS-68 (-2 to 40°C).
     * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
     * @return  Salinity PSS-78.
     *
     * @todo What to do with pression!?! Check the equation of state.
     */
    public static double salinity(double C, final double T, final double P) {
        C /= STANDARD_CONDUCTIVITY;
        if (!(C < 5E-4)) { // use '!' in order to accept NaN
            final double XR = sqrt(C/(poly(T,PSS78_C) * (1.0 + poly(P,PSS78_E) * P /
                              ((PSS78_D[1] * T+PSS78_D[0]) * T + 1.0 + (PSS78_D[3] * T + PSS78_D[2]) * C))));
            final double S  = sal(XR, T-15.0);  // Do not use an 'assert' statement invoking 'cond'.
            if (!(S>=42)) return S; // use '!' to accept NaN
            /*
             * Calcule la salinité pour une eau de conductivité,
             * de température et de pression données. Cet algorithme
             * doit être utilisé lorsque l'on s'attend à une salinité
             * entre 42 et 50.
             */
            return 35 * C + C * (C-1) * (poly(C,PSS78_AR) + T * (poly(T,PSS78_AT) + C *
                        (PSS78_CR[0] + PSS78_CR[1] * C + PSS78_CR[2] * T)));
            // TODO: VERIFIER CE QUE DEVIENT LA PRESSION ET IMPLEMENTER L'EQUATION D'ETAT.
        } else {
            return 0; // Zero conductivity trap
        }
    }

    /**
     * Computes conductivity as a function of salinity, temperature and pressure.
     *
     * @param S Salinity PSS-78 (0 to 42)
     * @param T Temperature ITS-68 (-2 to 40°C)
     * @param P Pressure (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
     * @return  Conductivity in mS/cm.
     */
    public static double conductivity(final double S, final double T, final double P) {
        if (!(S < 0.02)) { // use '!' in order to accept NaN
            double XT = T-15.0;
            double RT = sqrt(S / 35.0);  // First approximation
            double SI = sal(RT,XT);
            for (int n=0; n<10; n++) {    // Iteration loop begin here with a maximum of 10 cycles
                RT += (S-SI)/dsal(RT,XT);
                SI = sal(RT,XT);
                if (abs(SI-S) < 1E-4) {
                    break;
                }
            }
            double RTT = poly(T,PSS78_C)*(RT*RT);
            double AT  = PSS78_D[3]*T + PSS78_D[2];
            double BT  = (PSS78_D[1]*T + PSS78_D[0])*T + 1.0;
            double CP  = RTT*(BT + poly(P,PSS78_E)*P);
            BT -= RTT*AT;
            // Solve quadratic equation for C = RT35*RT*(1+C/AR+b)
            double cnd = 0.5*(sqrt(abs((BT*BT) + 4.0*AT*CP)) - BT)/AT;
            return cnd*STANDARD_CONDUCTIVITY;
        } else {
            return 0; // Zero salinity trap
        }
    }

    /**
     * Computes specific heat as a function of salinity, temperature and pressure.
     *
     * @param S Salinity PSS-78.
     * @param T Temperature (°C).
     * @param P Pressure (dbar), not including atmospheric pressure.
     * @return  Specific heat (J/(kg&times;°C)).
     */
    public static double specificHeat(final double S, final double T, double P) {
        P /= 10.0;
        final double SR = sqrt(S);
        return (poly(T,HEAT_CC) + (poly(T,HEAT_BB)*SR + poly(T,HEAT_AA))*S +
            (((poly(T,HEAT_C)*P + poly(T,HEAT_B) )*P + poly(T,HEAT_A) )*P) +
            ((((HEAT_J*SR+poly(T,HEAT_H))*S*P + (HEAT_G*SR+poly(T,HEAT_F))*S)*P +
                          (poly(T,HEAT_E)*SR+poly(T,HEAT_D))*S )*P));
    }

    /**
     * Computes fusion temperature (melting point) as a function of salinity and pressure.
     *
     * @param S Salinity PSS-78.
     * @param P Pressure (dbar), not including atmospheric pressure.
     * @return  Melting point (°C).
     */
    public static double fusionTemperature(final double S, final double P) {
        return (-0.0575 + 1.710523E-3*sqrt(S) + -2.154996E-4*S)*S + -7.53E-4*P;
    }

    /**
     * Computes adiabetic temperature gradient as a function of salinity, temperature and pressure.
     *
     * @param S Salinity PSS-78.
     * @param T Temperature (°C).
     * @param P Pressure (dbar), not including atmospheric pressure.
     * @return  Adiabetic temperature gradient (°C/dbar).
     */
    public static double adiabeticTemperatureGradient(double S, final double T, final double P) {
        S -= 35.0;
        return (poly(T,GRAD_A) + poly(T,GRAD_B)*S +
               (poly(T,GRAD_C) + poly(T,GRAD_D)*S + poly(T,GRAD_E)*P)*P);
    }

    /**
     * Computes depth as a function of pressure and latitude.
     *
     * @param  P Pressure (dbar), not including atmospheric pressure.
     * @param  lat Latitude in degrees (-90 to 90°)
     * @return Depth (m).
     */
    public static double depth(final double P, double lat) {
        lat = sin(lat);
        lat *= lat;
        lat = 9.780318*( 1.0 + 5.2788E-3*lat + 2.36E-5*(lat*lat));
        return poly(P,DEPTH_C)*P / (lat+(0.5*2.184E-6)*P);
    }

    /**
     * Computes sound velocity as a function of salinity, temperature and pressure.
     *
     * @param S Salinity PSS-78.
     * @param T Temperature (°C).
     * @param P Pressure (dbar), not including atmospheric pressure.
     * @return  Sound velocity (m/s).
     */
    public static double soundVelocity(final double S, final double T, final double P) {
        // S^0 terms
        final double CW = ((poly(T,SOUND_C3) *P + poly(T,SOUND_C2))*P +
                            poly(T,SOUND_C1))*P + poly(T,SOUND_C0);
        // S^1 terms
        final double A  = ((poly(T,SOUND_A3) *P + poly(T,SOUND_A2))*P +
                            poly(T,SOUND_A1))*P + poly(T,SOUND_A0);
        // S^3/2 terms
        final double B  = poly(T,SOUND_B0) + poly(T,SOUND_B1)*P;

        // S^2 terms
        final double D  = SOUND_D0 + SOUND_D1*P;

        // sound speed return
        return CW + (D*S + B*sqrt(S) + A)*S;
    }

    /**
     * Computes saturation in disolved oxygen as a function of salinity and temperature.
     *
     * @param S Salinity PSS-78.
     * @param T Temperature (°C).
     * @return  Saturation in disolved oxygen (µmol/kg).
     */
    public static double saturationO2(final double S, double T) {
        T += 273.15;
        return exp(poly_inv(T,O2_AT) + S*poly_inv(T,O2_AS));
    }

    /**
     * Calcule la valeur d'un polynôme.
     * Cette fonction calcule la valeur de:
     *
     * {@preformat java
     *    y = C[0] + C[1]*x + C[2]*x² + C[3]*x³
     * }
     *
     * où C est un vecteur de coéfficients transmis en argument.
     * Une exception sera levée si ce tableau ne contient pas
     * au moins 1 élément.
     *
     * @param x Valeur x à laquelle calculer le polynôme.
     * @param c Coéfficients C du polynôme.
     * @return  La valeur du polynôme au x spécifié.
     *
     * @see #poly_inv(double,double[])
     */
    private static double poly(final double x, final double[] c) {
        int n = c.length-1;
        double y = c[n];
        while (n > 0) {
            y = y*x + c[--n];
        }
        return y;
    }

    /**
     * Calcule la valeur de:
     *
     * {@preformat java
     *    y = C[0] + C[1]/x + C[2]/x² + C[3]/x³
     * }
     *
     * où C est un vecteur de coéfficients transmis en argument.
     * Une exception sera levée si ce tableau ne contient pas
     * au moins 1 élément.
     *
     * @param x Valeur x à laquelle calculer le polynôme.
     * @param C Coéfficients C du polynôme.
     * @return  La valeur du polynôme au x spécifié.
     *
     * @see #poly(double,double[])
     */
    private static double poly_inv(final double x, final double[] c) {
        int n = c.length-1;
        double y = c[n];
        while (n > 0) {
            y = y/x + c[--n];
        }
        return y;
    }
}
